function [EL_C2r,EA_C2r,EW_C2r,EL_C2,EL_T2,TE_2] = smom_function(param)

        %=========================================================================%
        % Part 1: parameters
        %=========================================================================%
        rng(1,'twister');
        s = rng;
 	    sigma = param(1);      %EIS      
	    beta  = .97;           %discount factor 
        aT    = 81;            %terminal age
        aS    = 50;            %start age
        r     = .029;          %interest rate
        pi    = .028;          %inflation
        V_term = 0;            %terminal value
        z_age_min = 56;        %earliest freeze age
        z_age_max = 64;        %latest freeze age
        N_obs  = 5000;         %Number of simulated individuals
        ages = aS:1:aT-1;      %ages
        t = (1:(aT-aS))';      %time periods
        scaler = 1;            %scaler on consumption utility 
        solver_param = [sigma;beta;aT;aS;r;V_term;scaler;pi]; %collect some params
      
        %-------------------------------------------------------------------------%
        %Calculate annuity value and earnings per period from HRS
        %-------------------------------------------------------------------------%
        %1. HRS db wealth by age
        pdata = csvread('compwealth_frz_w10_22.csv',1,0);
        %col2 is age
        %col3 is avg db wealth in 2010 dollars
        %col4 is avg earnings in 2010 dollars
        dbw = pdata(aS-48+1:aT-48,2:3);
        y_raw = pdata(aS-48+1:end,4);
        %Smooth the earnings profile using cubic polynomic fit
        y_fit_param =  polyfit(ages',y_raw,3);
        y = polyval(y_fit_param,ages');
        %percent earnings losses from freezes (periods 0-4)
        y_cut = [.037;.020;.015;.048;.016];
        %2. SS annuity by retirement age (hh-level)
        b_ss = csvread('ssw10.csv',1,0);
        b_ss = b_ss(aS-48+1:end,4);
        %3. SSA mortality rates by age (48+)
        pdeath = csvread('SSA_mortality_2010.csv',1,0);
        pdeath = pdeath(aS+1:end,:);
        avpdeath = mean(pdeath(:,2:3),2);
        %4. annuity factors
        afactor = zeros(length(avpdeath),1); 
        %loop over age
        for a = 1:length(afactor)
            disc = zeros(length(afactor)-a+1,1);
            disc(1) = 1;
            for i = 2:length(disc)
                disc(i) = (1-avpdeath(a+i-1))*(1/((1+r)*(1+pi)))^(i-1); 
            end

            afactor(a) = sum(disc);
        end     
        afactor = [pdeath(1:aT-aS,1) afactor(1:aT-aS)];
        %DB annuity by retirement age
        b_db = dbw(:,2)./afactor(:,2);
        %total annuity benefit
        b = b_db + b_ss;
        
        %-------------------------------------------------------------------------%
        % Normalized mortality rate
        %-------------------------------------------------------------------------%
        mortp = avpdeath(1:aT-aS);
        mortp_80 = mortp./avpdeath(aT-aS+1);
        
        %-------------------------------------------------------------------------%
        % Asset grids (A = non-pension; W = DC pension)
        %-------------------------------------------------------------------------%
        %A
        A_gridN = 20;
        split = 5;
        A_mins = [0;50000;150000;500000;1000000;2000000];
        As=zeros(A_gridN/split,split);
        for i = 1:5
            step = (A_mins(i+1)-A_mins(i))/(A_gridN/split);
            As(:,i) = A_mins(i):step:(A_mins(i+1)-step);
        end
        A = reshape(As,[A_gridN,1]);
        
        %W
        W_gridN = 20;
        split = 5;
        W_mins = [0;25000;75000;250000;500000;1500000];
        Ws=zeros(W_gridN/split,split);
        for i = 1:5
            step = (W_mins(i+1)-W_mins(i))/(W_gridN/split);
            Ws(:,i) = W_mins(i):step:(W_mins(i+1)-step);
        end
        W = reshape(Ws,[W_gridN,1]);    
        
        %-------------------------------------------------------------------------%
        % Average tax rates from nber taxsim (age-dependent grids)
        %-------------------------------------------------------------------------%       
        worktax = csvread('workinc_taxsim_a_20.csv',0,0);
        %col1 = age, col2 = earn_before_dc, col3 =fica, col4=earn_after_dc
        %col5=fiitax, col6=dc_contrib
        itau_w = worktax(:,5)./worktax(:,4); %income tax
        ptau_w = 0.5*worktax(:,3)./worktax(:,2); %worker share of FICA
      
        rettax = csvread('retinc_taxsim_a_20.csv',0,0);
        %col1 = fiitax, col2 = age, col3=pension distribution, col4 = SS
        itau_r = rettax(:,1)./(rettax(:,3)+rettax(:,4));
             
        %Collect parameters for model solution
        taxf = cell(3,1);
        taxf{1,1}=itau_w;
        taxf{2,1}=ptau_w;
        taxf{3,1}=itau_r;
        
        %freeze-affected average tax rates 
        taxfz = cell(9,3);
        zindex = 1;
        for aa = 56:64         
            i_z = find(ages==aa);
            
            %While working
            fname = strcat('workinc_taxsim_a_20_z',num2str(i_z),'.csv');
            worktax_z = csvread(fname,0,0);
            itau_w_z = worktax_z(:,5)./worktax_z(:,4); %income tax
            ptau_w_z = 0.5*worktax_z(:,3)./worktax_z(:,2); %worker share of FICA
            taxfz{zindex,1} = itau_w_z;
            taxfz{zindex,2} = ptau_w_z;
     
            %Retirement
            fname = strcat('retinc_taxsim_a_20_z',num2str(i_z),'.csv');
            rettax_z = csvread(fname,0,0);
            itau_r_z = rettax_z(:,1)./(rettax_z(:,3)+rettax_z(:,4));
            taxfz{zindex,3} = itau_r_z;
            zindex = zindex+1;
        end
        
        %-------------------------------------------------------------------------%
        %Disutility of work process 
        %-------------------------------------------------------------------------%
    	mu      = 0;  %mean normal underlying f
        rho     = 1/(1+exp(param(2)));  %persistence of random process 
    	sigma_v = exp(param(3)); %sd of normal underlying f
    	phi     = param(4); %age slope
        gamma   = param(5); %constant
    	f_N     = 6; %grid size

        %=========================================================================%
        % Part 2: solve model for value and policy function 
        % without any freeze activity
        %=========================================================================%

        %Cell array to store value and policy functions for each age
        optfunc = cell(length(t),8);

        %-------------------------------------------------------------------------%
        %Solve age T-1 policy and value functions outside the loop
        %-------------------------------------------------------------------------%	
        %Set last period=true
        last=1;
        a=0;
        [V_w,A_polw,W_polw]=vwork_5(A,W,last,a,y,optfunc,solver_param,...
                    f_N,mu,rho,sigma_v,phi,gamma,taxf,mortp_80);
        [V_r,A_polr,W_polr]=vret_4(A,W,last,a,b_db,b_ss,optfunc,solver_param,taxf,mortp_80);
	
        
        % Outer value function and labor supply choice
        V = zeros(size(V_w)); 
        L = zeros(size(V_w)); 
       
        for i = 1:length(A)
            for j = 1:length(W)
                for k = 1:f_N
                    V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                    L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                end
            end
        end

        %Store value and policy functions
        optfunc(1,:) = {V_w;A_polw;W_polw;V_r;A_polr;W_polr;V;L};
        
        %-------------------------------------------------------------------------%
        %Solve policy and value functions for remaning periods
        %-------------------------------------------------------------------------%      
        %Set last period = false
        last=0;
        
        %Loop through ages
        for a = 1:length(t)-1 
    

               [V_w,A_polw,W_polw]=vwork_5(A,W,last,a,y,optfunc,solver_param,...
                    f_N,mu,rho,sigma_v,phi,gamma,taxf,mortp_80);
               [V_r,A_polr,W_polr]=vret_4(A,W,last,a,b_db,b_ss,optfunc,solver_param,taxf,mortp_80);


		% Outer value function and labor supply choice
		V = zeros(size(V_w)); 
		L = zeros(size(V_w)); 

		for i = 1:length(A)
		    for j = 1:length(W)
			for k = 1:f_N
			    V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
			    L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
			end
		    end
		end

		%Store value and policy functions
        optfunc(a+1,:) = {V_w;A_polw;W_polw;V_r;A_polr;W_polr;V;L};

        end
        

        %=========================================================================%
        % Part 3: solve model for value and policy function 
        % separately for each potential freeze age 
        %=========================================================================%	
        %Cell array to store value and policy functions for each freeze age
        % dimension is age x freeze ages x number of functions (8)
        optfunc_frz = cell(length(t),z_age_max-z_age_min+1,8);
        lastfuncW0 = [];
        lastfuncR0 = [];
        
        parfor z = 1:z_age_max-z_age_min+1 %%loop over possible freeze ages
            
            %Freeze the path of DB pension annuities (keep SS wealth fixed)
            i_z = find(ages==z_age_min-1+z);
            b_db_z = b_db;
            b_db_z(i_z-1:end) = b_db(i_z-1); %no increment to db annuity at freeze age
            %freeze fixes nominal annuity; compute real annuity at each
            %prospective retirement age
            for i = i_z:length(b_db_z) 
                b_db_z(i) = b_db_z(i)*(1/(1+pi))^(i-i_z);
            end
            b_z = b_db_z + b_ss;
            %earn profile given freeze at age z
            y_loss = ones(length(y),1);
            y_loss(i_z:i_z+length(y_cut)-1) = 1-y_cut;
            y_z = y.*y_loss;

            %---------------------------------------------------------------------%
            %Solve age T-1 policy and value functions outside the age loop
            %---------------------------------------------------------------------%
            %Set last period=true
            last=1;
            a=0;

            [V_w,A_polw,W_polw]=vworkz_5(A,W,last,a,y_z,lastfuncW0,solver_param,...
                        f_N,mu,rho,sigma_v,phi,gamma,taxfz,mortp_80,z);
            [V_r,A_polr,W_polr]=vretz_4(A,W,last,a,b_db_z,b_ss,lastfuncR0,solver_param,taxfz,mortp_80,z);


            % Outer value function and labor supply choice
            V = zeros(size(V_w)); 
            L = zeros(size(V_w)); 

            for i = 1:length(A)
                for j = 1:length(W)
                    for k = 1:f_N
                        V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                        L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                    end
                end
            end

            %Store value and policy functions
            funcs = {V_w,A_polw,W_polw,V_r,A_polr,W_polr,V,L};   

            %-------------------------------------------------------------------------%
            %Solve policy and value functions for remaning periods
            %-------------------------------------------------------------------------%
            %Set last period = false
            last=0;

            %Loop through ages
            for a = 1:length(t)-1 

                lastfuncW = cell2mat(funcs(a,7));
                lastfuncR = funcs(a,4);

                [V_w,A_polw,W_polw]=vworkz_5(A,W,last,a,y_z,lastfuncW,solver_param,...
                        f_N,mu,rho,sigma_v,phi,gamma,taxfz,mortp_80,z);
                [V_r,A_polr,W_polr]=vretz_4(A,W,last,a,b_db_z,b_ss,lastfuncR,solver_param,taxfz,mortp_80,z);


                % Outer value function and labor supply choice
                V = zeros(size(V_w)); 
                L = zeros(size(V_w)); 

                for i = 1:length(A)
                    for j = 1:length(W)
                    for k = 1:f_N
                        V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                        L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                    end
                    end
                end

                %Store value and policy functions
                funcs0 = {V_w,A_polw,W_polw,V_r,A_polr,W_polr,V,L};
                funcs = [funcs;funcs0];

            end

            optfunc_frz(:,z,:) = funcs;

        end
        
       
        %=========================================================================%
        %Part 4: Simulate shocks and solve for labor supply and saving
        %=========================================================================%
        
        %Cell arrays to hold simulated paths 
        LsimC = cell(z_age_max-z_age_min+1,1);
        AsimC = cell(z_age_max-z_age_min+1,1);
        WsimC = cell(z_age_max-z_age_min+1,1);
        LsimT = cell(z_age_max-z_age_min+1,1);
        AsimT = cell(z_age_max-z_age_min+1,1);
        WsimT = cell(z_age_max-z_age_min+1,1);
      
        %Nodes of discrete AR(1) process -- need max and min of these to
        %assign values out of bounds in the simulation
        [f,~] = rouwenhorst(f_N,mu,rho,sigma_v);
        hor = 12; %horizon
        
        %Simulated observations per age (prop to LFP in HRS)
        age_shares = csvread('age_shares.csv',1,1);
        obs = round(N_obs.*age_shares(:,1));
        
        %Disutility shock process
        D = cell(9,1);
        parfor a = 1:9
            %simulate AR(1) shocks. 1000 period burn-in... 
            f_sim0 = zeros(1000+(hor+6),obs(a));
            rng(s);
            v = normrnd(0,sigma_v,[length(f_sim0),obs(a)]);
            for i = 2:length(f_sim0)
                f_sim0(i,:) = (1-rho)*mu+rho*f_sim0(i-1,:)+v(i,:);
            end
            %Trim away burn-in periods
            f_sim0 = f_sim0(end-(hor+6-1):end,:);
            D{a,1} = f_sim0;
        end
       
        %Initial assets --moments from HRS
        A0 = cell(9,1);
        W0 = cell(9,1);
        A_params = csvread('A_params.csv',1,0);
        J_params = csvread('J_params.csv',1,0);
        for a = 1:9
            rng(s);
            %Number of individuals with 0 DC wealth
            frac0obs = round(A_params(a,4)*obs(a)); 
            %Non-DC wealth for individuals with 0 DC wealth
            A_temp0 = lognrnd(A_params(a,2),A_params(a,3),[frac0obs,1]);
            %DC wealth for individuals with no DC balances = 0
            W_temp0 = zeros(frac0obs,1);
            %Joint wealth distribution for individuals with >0 DC wealth
            mvnmu = J_params(a,2:3);
            mvnsig = [(J_params(a,4))^2 J_params(a,6);J_params(a,6) (J_params(a,5))^2];
            J_wealth = exp(mvnrnd(mvnmu,mvnsig,(obs(a)-frac0obs)));
            W_temp1 = J_wealth(:,1);
            A_temp1 = J_wealth(:,2);
            A0{a,1} = [A_temp0' A_temp1'];
            W0{a,1} = [W_temp0' W_temp1'];
        end

        %use policy functions to obtain optimal labor supply and asset
        %accumulation for each simulated individual in each period 
        %-------------------------------------------------------------------------%
        %Control group
        %-------------------------------------------------------------------------%
        
        parfor p = 1:9 
            L_sim = zeros(hor+6,obs(p));
            A_sim = zeros(hor+6,obs(p));
            W_sim = zeros(hor+6,obs(p));
            A_sim(1,:) = A0{p,1};
            W_sim(1,:) = W0{p,1};
        
            %Set disutility distribution
            f_sim = D{p,1};
 
            %Employment choice for t = -5 outside the loop (t = period
            %relative to freeze) a = 1 when t = -5
            for i = 1:obs(p)
                a=1;    
                Afunc_w = optfunc{(aT-aS)-(p+a-1),2};
                Wfunc_w = optfunc{(aT-aS)-(p+a-1),3};      
                Afunc_r1 = optfunc{(aT-aS)-(p+a-1),5};
                Afunc_r  = Afunc_r1{p+a};
                Wfunc_r1 = optfunc{(aT-aS)-(p+a-1),6};  
                Wfunc_r  = Wfunc_r1{p+a};
                Lfunc = optfunc{(aT-aS)-(p+a-1),8};
                %Put extreme draws of A,W,f onto the grid 
                %A
                if A_sim(a,i)>max(A)
                    A_sim(a,i) = max(A);
                end
                %W
                if W_sim(a,i)>max(W)
                    W_sim(a,i) = max(W);
                end	
                %f                
                if f_sim(a,i)<min(f) 
                    f_sim(a,i) = min(f);
                elseif f_sim(a,i)>max(f)
                    f_sim(a,i) = max(f);
                end      

                %Simulated values 
                L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Lfunc,A_sim(a,i),W_sim(a,i),f_sim(a,i)));

                if L_sim(a,i) == 1
                    A_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Afunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                    W_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Wfunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                elseif L_sim(a,i) == 0 
                    A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                    W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));   
                end  
            end

            for a = 2:hor+5 
               for i = 1:obs(p)

                    Afunc_w = optfunc{(aT-aS)-(p+a-1),2};
                    Wfunc_w = optfunc{(aT-aS)-(p+a-1),3};      
                    Afunc_r1 = optfunc{(aT-aS)-(p+a-1),5};
                    Afunc_r  = Afunc_r1{p+a};
                    Wfunc_r1 = optfunc{(aT-aS)-(p+a-1),6};  
                    Wfunc_r  = Wfunc_r1{p+a};                  
                    Lfunc = optfunc{(aT-aS)-(p+a-1),8};

                    %Put extreme draws of A,W,f onto the grid 
                    %A
                    if A_sim(a,i)>max(A)
                        A_sim(a,i) = max(A);
                    end
                    %W
                    if W_sim(a,i)>max(W)
                        W_sim(a,i) = max(W);
                    end	
                    %f                
                    if f_sim(a,i)<min(f) 
                        f_sim(a,i) = min(f);
                    elseif f_sim(a,i)>max(f)
                        f_sim(a,i) = max(f);
                    end        

                    %Enforce self-absorbing retirement (value functions impose it)
                    if  L_sim(a-1,i) == 1                 
                        L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Lfunc,A_sim(a,i),W_sim(a,i),f_sim(a,i)));                   
                            if L_sim(a,i) == 1
                                A_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Afunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                                W_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Wfunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                            elseif L_sim(a,i) == 0 
                                A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                                W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));   
                            end                    
                    elseif L_sim(a-1,i) == 0
                        L_sim(a,i) = 0;
                        A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                        W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));
                    end

                end
            end
            
            %Final period employment choice (A,W already determined)
            a = hor+6;
            for i = 1:obs(p)

                    Afunc_w = optfunc{(aT-aS)-(p+a-1),2};
                    Wfunc_w = optfunc{(aT-aS)-(p+a-1),3};      
                    Afunc_r1 = optfunc{(aT-aS)-(p+a-1),5};
                    Afunc_r  = Afunc_r1{p+a};
                    Wfunc_r1 = optfunc{(aT-aS)-(p+a-1),6};  
                    Wfunc_r  = Wfunc_r1{p+a};                  
                    Lfunc = optfunc{(aT-aS)-(p+a-1),8};

                    %Put extreme draws of A,W,f onto the grid 
                    %A
                    if A_sim(a,i)>max(A)
                        A_sim(a,i) = max(A);
                    end
                    %W
                    if W_sim(a,i)>max(W)
                        W_sim(a,i) = max(W);
                    end	
                    %f                
                    if f_sim(a,i)<min(f) 
                        f_sim(a,i) = min(f);
                    elseif f_sim(a,i)>max(f)
                        f_sim(a,i) = max(f);
                    end        

                    %Enforce self-absorbing retirement (value functions impose it)
                    if  L_sim(a-1,i) == 1                 
                        L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Lfunc,A_sim(a,i),W_sim(a,i),f_sim(a,i)));                                     
                    elseif L_sim(a-1,i) == 0
                        L_sim(a,i) = 0;
                    end

            end
        
            LsimC{p,1} = L_sim;
            AsimC{p,1} = A_sim;
            WsimC{p,1} = W_sim;
        end
    
    %-------------------------------------------------------------------------%
	%Treated (i.e. DB frozen) group
	%-------------------------------------------------------------------------%
    
     parfor p = 1:9 
        L_sim = zeros(hor+6,obs(p));
        A_sim = zeros(hor+6,obs(p));
        W_sim = zeros(hor+6,obs(p));
        A_sim(1,:) = A0{p,1};
        W_sim(1,:) = W0{p,1};

        %Set disutility distribution
        f_sim = D{p,1};

        %Employment choice for t = -5 outside the loop (t = period
        %relative to freeze) a = 1 when t = -5
        for i = 1:obs(p)
            a=1;    
	    	
            %No-freeze policy functions
	    	Afunc_w = optfunc{(aT-aS)-(p+a-1),2};
            Wfunc_w = optfunc{(aT-aS)-(p+a-1),3};      
            Afunc_r1 = optfunc{(aT-aS)-(p+a-1),5};
            Afunc_r  = Afunc_r1{p+a};
            Wfunc_r1 = optfunc{(aT-aS)-(p+a-1),6};  
            Wfunc_r  = Wfunc_r1{p+a};                  
            Lfunc = optfunc{(aT-aS)-(p+a-1),8};

            %Put extreme draws of A,W,f onto the grid 
            %A
            if A_sim(a,i)>max(A)
                A_sim(a,i) = max(A);
            end
            %W
            if W_sim(a,i)>max(W)
                W_sim(a,i) = max(W);
            end	
            %f                
            if f_sim(a,i)<min(f) 
                f_sim(a,i) = min(f);
            elseif f_sim(a,i)>max(f)
                f_sim(a,i) = max(f);
            end      

                %Simulated values 
                L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Lfunc,A_sim(a,i),W_sim(a,i),f_sim(a,i)));
                if L_sim(a,i) == 1
                    A_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Afunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                    W_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Wfunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                elseif L_sim(a,i) == 0 
                    A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                    W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));   
                end  
        end
        
        for a = 2:5 % pre-freeze after initial period
           for i = 1:obs(p)

                %No-freeze policy functions
                Afunc_w = optfunc{(aT-aS)-(p+a-1),2};
                Wfunc_w = optfunc{(aT-aS)-(p+a-1),3};      
                Afunc_r1 = optfunc{(aT-aS)-(p+a-1),5};
                Afunc_r  = Afunc_r1{p+a};
                Wfunc_r1 = optfunc{(aT-aS)-(p+a-1),6};  
                Wfunc_r  = Wfunc_r1{p+a};                 
                Lfunc = optfunc{(aT-aS)-(p+a-1),8};

                %Put extreme draws of A,W,f onto the grid 
                %A
                if A_sim(a,i)>max(A)
                    A_sim(a,i) = max(A);
                end
                %W
                if W_sim(a,i)>max(W)
                    W_sim(a,i) = max(W);
                end	
                %f                
                if f_sim(a,i)<min(f) 
                    f_sim(a,i) = min(f);
                elseif f_sim(a,i)>max(f)
                    f_sim(a,i) = max(f);
                end        

                %Enforce self-absorbing retirement (value functions impose it)
                if  L_sim(a-1,i) == 1                 
                    L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Lfunc,A_sim(a,i),W_sim(a,i),f_sim(a,i)));                   
                        if L_sim(a,i) == 1
                            A_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Afunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                            W_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Wfunc_w,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                        elseif L_sim(a,i) == 0 
                            A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                            W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));   
                        end                    
                elseif L_sim(a-1,i) == 0
                    L_sim(a,i) = 0;
                    A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_r,A_sim(a,i),W_sim(a,i));
                    W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_r,A_sim(a,i),W_sim(a,i));
                end
                    
            end
        end
        

        for a = 6:hor+5 %after freeze
           for i = 1:obs(p)

                %Freeze policy functions
                Afunc_wZ = optfunc_frz{(aT-aS)-(p+a-1),p,2};
                Wfunc_wZ = optfunc_frz{(aT-aS)-(p+a-1),p,3};      
                Afunc_rZ1 = optfunc_frz{(aT-aS)-(p+a-1),p,5};
                Afunc_rZ  = Afunc_rZ1{p+a};
                Wfunc_rZ1 = optfunc_frz{(aT-aS)-(p+a-1),p,6};  
                Wfunc_rZ  = Wfunc_rZ1{p+a};
                LfuncZ = optfunc_frz{(aT-aS)-(p+a-1),p,8};

                %Put extreme draws of A,W,f onto the grid 
                %A
                if A_sim(a,i)>max(A)
                    A_sim(a,i) = max(A);
                end
                %W
                if W_sim(a,i)>max(W)
                    W_sim(a,i) = max(W);
                end	
                %f                
                if f_sim(a,i)<min(f) 
                    f_sim(a,i) = min(f);
                elseif f_sim(a,i)>max(f)
                    f_sim(a,i) = max(f);
                end        

                %Enforce self-absorbing retirement (value functions impose it)
                if  L_sim(a-1,i) == 1                 
                    L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,LfuncZ,A_sim(a,i),W_sim(a,i),f_sim(a,i)));                   
                        if L_sim(a,i) == 1
                            A_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Afunc_wZ,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                            W_sim(a+1,i) = lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,Wfunc_wZ,A_sim(a,i),W_sim(a,i),f_sim(a,i));
                        elseif L_sim(a,i) == 0 
                            A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_rZ,A_sim(a,i),W_sim(a,i));
                            W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_rZ,A_sim(a,i),W_sim(a,i));   
                        end                    
                elseif L_sim(a-1,i) == 0
                    L_sim(a,i) = 0;
                    A_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Afunc_rZ,A_sim(a,i),W_sim(a,i));
                    W_sim(a+1,i) = lininterp2(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),Wfunc_rZ,A_sim(a,i),W_sim(a,i));
                end

            end
        end
        
        %Final period employment choice (A,W already determined)
        a = hor+6;
        for i = 1:obs(p)
            %Freeze policy functions
            Afunc_wZ = optfunc_frz{(aT-aS)-(p+a-1),p,2};
            Wfunc_wZ = optfunc_frz{(aT-aS)-(p+a-1),p,3};      
            Afunc_rZ1 = optfunc_frz{(aT-aS)-(p+a-1),p,5};
            Afunc_rZ  = Afunc_rZ1{p+a};
            Wfunc_rZ1 = optfunc_frz{(aT-aS)-(p+a-1),p,6};  
            Wfunc_rZ  = Wfunc_rZ1{p+a};
            LfuncZ = optfunc_frz{(aT-aS)-(p+a-1),p,8};

            %Put extreme draws of A,W,f onto the grid 
            %A
            if A_sim(a,i)>max(A)
                A_sim(a,i) = max(A);
            end
            %W
            if W_sim(a,i)>max(W)
                W_sim(a,i) = max(W);
            end	
            %f                
            if f_sim(a,i)<min(f) 
                f_sim(a,i) = min(f);
            elseif f_sim(a,i)>max(f)
                f_sim(a,i) = max(f);
            end        

            %Enforce self-absorbing retirement (value functions impose it)
            if  L_sim(a-1,i) == 1                 
                L_sim(a,i) = round(lininterp3(A.*(1+r)^(p+a-1),W.*(1+r)^(p+a-1),f,LfuncZ,A_sim(a,i),W_sim(a,i),f_sim(a,i)));                                  
            elseif L_sim(a-1,i) == 0
                L_sim(a,i) = 0;
            end
            
        end

        LsimT{p,1} = L_sim;
        AsimT{p,1} = A_sim;
        WsimT{p,1} = W_sim;
        
     end

    %-------------------------------------------------------------------------%
    % Compute moments by averaging over simulated individuals
    %-------------------------------------------------------------------------%

    %Concatenate control group       
    L_C2 = [LsimC{1,1} LsimC{2,1} LsimC{3,1} LsimC{4,1} LsimC{5,1} LsimC{6,1} ...
            LsimC{7,1} LsimC{8,1} LsimC{9,1}];
    A_C2 = [AsimC{1,1} AsimC{2,1} AsimC{3,1} AsimC{4,1} AsimC{5,1} AsimC{6,1} ...
            AsimC{7,1} AsimC{8,1} AsimC{9,1}];
    W_C2 = [WsimC{1,1} WsimC{2,1} WsimC{3,1} WsimC{4,1} WsimC{5,1} WsimC{6,1} ...
            WsimC{7,1} WsimC{8,1} WsimC{9,1}];
       
    %Condition on sample of workers who are still in LF 5 years pre-freeze
    indx_1 = find(L_C2(1,:)==0);

    L_C2_r1 = L_C2;
    L_C2_r1(:,indx_1)=[];
    A_C2_r1 = A_C2;
    A_C2_r1(:,indx_1)=[];
    W_C2_r1 = W_C2;
    W_C2_r1(:,indx_1)=[];

    %Extract treated group data        
    L_T2 = [LsimT{1,1} LsimT{2,1} LsimT{3,1} LsimT{4,1} LsimT{5,1} LsimT{6,1} ...
            LsimT{7,1} LsimT{8,1} LsimT{9,1}];
    A_T2 = [AsimT{1,1} AsimT{2,1} AsimT{3,1} AsimT{4,1} AsimT{5,1} AsimT{6,1} ...
            AsimT{7,1} AsimT{8,1} AsimT{9,1}];
    W_T2 = [WsimT{1,1} WsimT{2,1} WsimT{3,1} WsimT{4,1} WsimT{5,1} WsimT{6,1} ...
            WsimT{7,1} WsimT{8,1} WsimT{9,1}];  
    
    indx_2 = find(L_T2(1,:)==0);

    L_T2_r1 = L_T2;
    L_T2_r1(:,indx_2)=[];
    
    %Make simulated data moments
    %ensure that the data matrix is not empty
    %1. Employment rates conditional on working 5 years pre-freeze
    if isempty(L_C2_r1)==0
        EL_C2r = mean(L_C2_r1(2:end,:),2); %mean employment rates
        EA_C2r = mean(A_C2_r1(2:end,:),2); %mean non-pension assets
        EW_C2r = mean(W_C2_r1(2:end,:),2); %mean DC wealth
    else
        EL_C2r = ones(hor+5,1); 
        EA_C2r = zeros(hor+5,1); 
        EW_C2r = zeros(hor+5,1); 
    end    
    %2. employment rates and ATE of freeze conditional on working 5 yrs pre freeze
    %age
    if isempty(L_C2_r1)==0 && isempty(L_T2_r1)==0
        EL_C2 = mean(L_C2_r1(6:end,:),2); %mean ER control
        EL_T2 = mean(L_T2_r1(6:end,:),2); %mean ER affected
        TE_2= (EL_T2-EL_C2); %Teffect (pp  diff in ER's)
    else
        EL_C2 = ones(hor+5,1); 
        EL_T2 = zeros(hor+5,1); 
        TE_2= zeros(hor+1,1);
    end
        
end      

